      subroutine solsanm

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'setls.h'
      include 'sanm.h'
      include 'eigv.h'
      include 'sanmjin.h'
      
      integer det
      real(NBF) cp011,cp012,cp021,cp022,cp111,cp112,cp121,cp122,cp211,cp212 &
     ,         cp221,cp222,cp311,cp312,cp321,cp322,cp411,cp412,cp421,cp422 &
     ,         krt11,krt12,krt21,krt22,rkrt11,rkrt12,rkrt21,rkrt22 &
     ,         ksq11,ksq12,ksq21,ksq22,rksq11,rksq12,rksq21,rksq22 &
     ,         skrt11,skrt12,skrt21,skrt22,ckrt11,ckrt12,ckrt21,ckrt22 &
     ,         eig11,eig12,eig21,eig22,beta11,beta12,beta21,beta22 &
     ,         rskrt11,rskrt12,rskrt21,rskrt22,a0,b0,c0,d0,a1,b1,c1,d1 &
     ,         rdet1,rdet2,const1,const2,const3,const4,const5,const6,const7,const8 &
     ,         temp1,temp2,temp3,temp4,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12 &
     ,         low,up11,up12,up21,up22,error1,error2,temp
!      real(NBF),dimension(:) :: j1_err(ng),j2_err(ng),j1_old(ng),j2_old(ng) 

!#define test
#ifdef test
      qs(0,1,1)=0.321510000
      qs(1,1,1)=0
      qs(2,1,1)=-0.002242270
      qs(3,1,1)=0
      qs(4,1,1)=-0.103879624
      qs(0,2,1)=0.260570000
      qs(1,2,1)=-0.113954400
      qs(2,2,1)=-0.129615633
      qs(3,2,1)=0.006926740
      qs(4,2,1)=0.266604857
      avgflx(1,1,1)=10.38437
      avgflx(1,2,1)=7.718542
      avgflx(2,1,1)=-0.00706014
      avgflx(2,2,1)=-0.0375938

      qs(0,1,2)=0.215290000
      qs(1,1,2)=0
      qs(2,1,2)=-0.016590801
      qs(3,1,2)=0
      qs(4,1,2)=0.011953675
      qs(0,2,2)=0.12666
      qs(1,2,2)=-0.0563346
      qs(2,2,2)=-0.007480921
      qs(3,2,2)=0.001116356
      qs(4,2,2)=-0.00863688
      avgflx(1,1,2)=2.1103759
      avgflx(1,2,2)=0.5986204
      avgflx(2,1,2)=-0.0058428
      avgflx(2,2,2)=0.0057008
#endif

! Update the source coefficients with the average source from CMFD
#define srcup
#ifdef srcup
      do im=1,ng
        do i=1,tnode
          do k=1,4
            if (qs(0,i,im).ne.0) then
              qs(k,i,im)=src(1,i,im)*(qs(k,i,im)/qs(0,i,im))
            endif              
          enddo        
          qs(0,i,im)=src(1,i,im)
        enddo
      enddo
#endif
      if (bcb.ne.0.or.bcu.ne.0) then
        det=1
      else
        det=0
      endif

!      do k=1,2
!        do i=1,tnode-1,2
        do i=1,tnode-1
          j=i+1
          j1_old=0
          j2_old=0
!
          if (det.ne.0.and.i.eq.1) then
            call sanm1n(i,j)
          endif

          do iin=1,100
            fsc=0
            do im=1,ng
! substitutions
            krt11=krt(1,i,im)
            krt12=krt(1,j,im)
            krt21=krt(2,i,im)
            krt22=krt(2,j,im)
            ksq11=ksq(1,i,im)
            ksq12=ksq(1,j,im)
            ksq21=ksq(2,i,im)
            ksq22=ksq(2,j,im)
            rkrt11=rkrt(1,i,im)
            rkrt12=rkrt(1,j,im)
            rkrt21=rkrt(2,i,im)
            rkrt22=rkrt(2,j,im)
            rksq11=rksq(1,i,im)
            rksq12=rksq(1,j,im)
            rksq21=rksq(2,i,im)
            rksq22=rksq(2,j,im)

            skrt11=sinh(krt11)
            skrt12=sinh(krt12)
            skrt21=sinh(krt21)
            skrt22=sinh(krt22)
            ckrt11=cosh(krt11)
            ckrt12=cosh(krt12)
            ckrt21=cosh(krt21)
            ckrt22=cosh(krt22)

            eig11=eigvec1(i,im)
            eig12=eigvec1(j,im)
            eig21=eigvec2(i,im)
            eig22=eigvec2(j,im)
            beta11=betax(1,i,im)
            beta12=betax(1,j,im)
            beta21=betax(2,i,im)
            beta22=betax(2,j,im)

            rskrt11=1/sinh(krt11)
            rskrt12=1/sinh(krt12)
            rskrt21=1/sinh(krt21)
            rskrt22=1/sinh(krt22)
! the coefficients of the particular solution
! cp/order/moment/l-rnode/
            cp011=imq1(i,im)*(5*qs(4,i,im)*(21+2*ksq(1,i,im))+qs(0,i,im)*kqu(1,i,im) &
                 +3*qs(2,i,im)*ksq(1,i,im))*rksx(1,i,im)
            cp012=imq1(j,im)*(5*qs(4,j,im)*(21+2*ksq(1,j,im))+qs(0,j,im)*kqu(1,j,im) &
                 +3*qs(2,j,im)*ksq(1,j,im))*rksx(1,j,im)
            cp021=imq2(i,im)*(5*qs(4,i,im)*(21+2*ksq(2,i,im))+qs(0,i,im)*kqu(2,i,im) &
                 +3*qs(2,i,im)*ksq(2,i,im))*rksx(2,i,im)
            cp022=imq2(j,im)*(5*qs(4,j,im)*(21+2*ksq(2,j,im))+qs(0,j,im)*kqu(2,j,im) &
                 +3*qs(2,j,im)*ksq(2,j,im))*rksx(2,j,im)
            cp111=imq1(i,im)*(qs(1,i,im)*ksq(1,i,im)+15*qs(3,i,im))*rkqu(1,i,im)
            cp112=imq1(j,im)*(qs(1,j,im)*ksq(1,j,im)+15*qs(3,j,im))*rkqu(1,j,im)
            cp121=imq2(i,im)*(qs(1,i,im)*ksq(2,i,im)+15*qs(3,i,im))*rkqu(2,i,im)
            cp122=imq2(j,im)*(qs(1,j,im)*ksq(2,j,im)+15*qs(3,j,im))*rkqu(2,j,im)
            cp211=imq1(i,im)*(qs(2,i,im)*ksq(1,i,im)+35*qs(4,i,im))*rkqu(1,i,im)
            cp212=imq1(j,im)*(qs(2,j,im)*ksq(1,j,im)+35*qs(4,j,im))*rkqu(1,j,im)
            cp221=imq2(i,im)*(qs(2,i,im)*ksq(2,i,im)+35*qs(4,i,im))*rkqu(2,i,im)
            cp222=imq2(j,im)*(qs(2,j,im)*ksq(2,j,im)+35*qs(4,j,im))*rkqu(2,j,im)
            cp311=imq1(i,im)*qs(3,i,im)*rksq(1,i,im)
            cp312=imq1(j,im)*qs(3,j,im)*rksq(1,j,im)
            cp321=imq2(i,im)*qs(3,i,im)*rksq(2,i,im)
            cp322=imq2(j,im)*qs(3,j,im)*rksq(2,j,im)
            cp411=imq1(i,im)*qs(4,i,im)*rksq(1,i,im)
            cp412=imq1(j,im)*qs(4,j,im)*rksq(1,j,im)
            cp421=imq2(i,im)*qs(4,i,im)*rksq(2,i,im)
            cp422=imq2(j,im)*qs(4,j,im)*rksq(2,j,im)

! Coefficients of the Cosh term
            cosal(1)=-(rskrt11*krt11*(eig11*cp011-eig21*cp011 &
                    -avgflx(1,i,im)+eig21*avgflx(2,i,im)))/(eig11-eig21)      
            cosar(1)=-(rskrt12*krt12*(eig12*cp012-eig22 &
                    *cp012-avgflx(1,j,im)+eig22*avgflx(2,j,im)))/(eig12-eig22)      
            cosal(2)=(rskrt21*krt21*(eig21*cp021-avgflx(1,i,im) &
                    +eig11*(-cp021+avgflx(2,i,im))))/(eig11-eig21)
            cosar(2)=(rskrt22*krt22*(eig22*cp022-avgflx(1,j,im) &
                    +eig12*(-cp022+avgflx(2,j,im))))/(eig12-eig22)
! substitutions
            a0=eig11*skrt11
            b0=eig12*skrt12
            c0=-ckrt11*krt11*beta11*(2+eig11)
            d0=ckrt12*krt12*beta12*(2+eig12)
            rdet1=1/(b0*c0-a0*d0)
            a1=skrt21
            b1=skrt22
            c1=-ckrt21*krt21*beta21
            d1=ckrt22*krt22*beta22
            rdet2=1/(b1*c1-a1*d1)

            const1=eig21*skrt21
            const2=eig22*skrt22
            const3=ckrt21*krt21*beta11*(2+eig21)
            const4=ckrt22*krt22*beta12*(2+eig22)
            const5=skrt11
            const6=skrt12
            const7=ckrt11*krt11*beta21
            const8=ckrt12*krt12*beta22

! fix
            ccpl11=cp011-cp111+cp211-cp311+cp411
            ccpl21=cp021-cp121+cp221-cp321+cp421
            ccpl12=cp011+cp111+cp211+cp311+cp411 
            ccpl22=cp021+cp121+cp221+cp321+cp421
            ccpr11=cp012-cp112+cp212-cp312+cp412 
            ccpr21=cp022-cp122+cp222-cp322+cp422
            ccpr12=cp012+cp112+cp212+cp312+cp412
            ccpr22=cp022+cp122+cp222+cp322+cp422
            ccpt11=cp111+3*cp211+6*cp311+10*cp411
            ccpt21=cp121+3*cp221+6*cp321+10*cp421
            ccpt12=cp112-3*cp212+6*cp312-10*cp412
            ccpt22=cp122-3*cp222+6*cp322-10*cp422

            temp1=eig11*(-cosal(1)*ckrt11-ccpl12)+eig12*(cosar(1)*ckrt12+ccpr11) &
                 +eig21*(-cosal(2)*ckrt21-ccpl22)+eig22*(cosar(2)*ckrt22+ccpr21)
           temp2=beta11*((eig11+2)*(cosal(1)*krt11*skrt11+ccpt11) &
                 +(eig21+2)*(cosal(2)*krt21*skrt21+ccpt21)) &
                 +beta12*((eig12+2)*(cosar(1)*krt12*skrt12-ccpt12) &
                 +(eig22+2)*(cosar(2)*krt22*skrt22-ccpt22))
           temp3=-cosal(1)*ckrt11+cosar(1)*ckrt12-cosal(2)*ckrt21 &
                 +cosar(2)*ckrt22-ccpl12+ccpr11-ccpl22+ccpr21
           temp4=beta21*(cosal(1)*krt11*skrt11+cosal(2)*krt21*skrt21+ccpt11+ccpt21) &
                 +beta22*(cosar(1)*krt12*skrt12+cosar(2)*krt22*skrt22-ccpt12-ccpt22)

            k1=(b0*const3+d0*const1)*rdet1
            k2=(-b0*const4+d0*const2)*rdet1
            k3=(-d0*temp1+b0*temp2)*rdet1
            k4=(c0*const1+a0*const3)*rdet1
            k5=(c0*const2-a0*const4)*rdet1
            k6=(-c0*temp1+a0*temp2)*rdet1
            k7=(b1*const7+const5*d1)*rdet2
            k8=(-b1*const8+const6*d1)*rdet2
            k9=(-d1*temp3+b1*temp4)*rdet2
            k10=(c1*const5+a1*const7)*rdet2
            k11=(c1*const6-a1*const8)*rdet2
            k12=(-c1*temp3+a1*temp4)*rdet2

            low=1/(1-k1*k7-k11*(k5+k2*k4*k7-k1*k5*k7)+k4*k8+k10*(k2+k2*k4*k8-k1*k5*k8))
            up11=k3-k11*k3*k5+k11*k2*k6+k3*k4*k8-k1*k6*k8-k12*(k2+k2*k4*k8-k1*k5*k8) &
                +k1*k9+k11*k2*k4*k9-k1*k11*k5*k9
            up12=-(k6+k3*k4*k7-k1*k6*k7-k12*(k5+k2*k4*k7-k1*k5*k7)+k4*k9 &
                +k10*(-k3*k5+k2*k6+k2*k4*k9-k1*k5*k9))
            up21=k11*k2*k6*k7-k6*k8-k10*k2*k6*k8+k12*(-k2*k7+k5*k8) &
                +k3*(k7-k11*k5*k7+k10*k5*k8)+k9+k10*k2*k9-k11*k5*k9
            up22=k12*(-1+k1*k7-k4*k8)-k10*(k3+k3*k4*k8-k1*k6*k8+k1*k9) &
               +k11*(k6+k3*k4*k7-k1*k6*k7+k4*k9)

! Coefficients of the Sinh term
            sinbl(1)=up11*low  
            sinbr(1)=up12*low   
            sinbl(2)=up21*low  
            sinbr(2)=up22*low

! Temporary coefficients storage
#define tempc
#ifdef tempc
            tcp01(i,im)=cp011
            tcp01(j,im)=cp012
            tcp11(i,im)=cp111
            tcp11(j,im)=cp112
            tcp21(i,im)=cp211
            tcp21(j,im)=cp212
            tcp31(i,im)=cp311
            tcp31(j,im)=cp312
            tcp41(i,im)=cp411
            tcp41(j,im)=cp412
            tcp02(i,im)=cp021
            tcp02(j,im)=cp022
            tcp12(i,im)=cp121
            tcp12(j,im)=cp122
            tcp22(i,im)=cp221
            tcp22(j,im)=cp222
            tcp32(i,im)=cp321
            tcp32(j,im)=cp322
            tcp42(i,im)=cp421
            tcp42(j,im)=cp422

            tcosa1(i,im)=cosal(1)
            tcosa1(j,im)=cosar(1)
            tcosa2(i,im)=cosal(2)
            tcosa2(j,im)=cosar(2)

            tsinb1(i,im)=sinbl(1)
            tsinb1(j,im)=sinbr(1)
            tsinb2(i,im)=sinbl(2)
            tsinb2(j,im)=sinbr(2)
#endif

! Surface flux
            if (i.eq.1) then
              sflux(1,i,im)=eig11*(ckrt11*cosal(1)-skrt11*sinbl(1)+ccpl11) &
                           +eig21*(ckrt21*cosal(2)-skrt21*sinbl(2)+ccpl21)
             sflux(2,i,im)=ckrt11*cosal(1)+ckrt21*cosal(2)-skrt11*sinbl(1) &
                           -skrt21*sinbl(2)+ccpl11+ccpl21
           elseif (j.eq.tnode) then
             sflux(1,j+1,im)=eig12*(ckrt12*cosar(1)+skrt12*sinbr(1)+ccpr12) &
                             +eig22*(ckrt22*cosar(2)+skrt22*sinbr(2)+ccpr22)
             sflux(2,j+1,im)=ckrt12*cosar(1)+ckrt22*cosar(2)+skrt12*sinbr(1) &
                             +skrt22*sinbr(2)+ccpr12+ccpr22
           endif
           sflux(1,j,im)=eig12*(-skrt12*sinbr(1)+ckrt12*cosar(1)+ccpr11) &
                         +eig22*(-skrt22*sinbr(2)+ckrt22*cosar(2)+ccpr21)
           sflux(2,j,im)=ckrt12*cosar(1)+ckrt22*cosar(2)-skrt12*sinbr(1)-skrt22*sinbr(2) &
                         +ccpr11+ccpr21

! Update the surface current
            jnet(1,j,im)=-beta12*((2+eig12)*krt12*(sinbr(1)*ckrt12-cosar(1)*skrt12) &
                        +(2+eig22)*krt22*(sinbr(2)*ckrt22-cosar(2)*skrt22) &
                        +eig12*ccpt12+eig22*ccpt22+2*(ccpt12+ccpt22)) 
            jnet(2,j,im)=-beta22*(sinbr(1)*ckrt12*krt12+sinbr(2)*ckrt22*krt22 &
                        -cosar(1)*skrt12*krt12-cosar(2)*skrt22*krt22+ccpt12+ccpt22)

! 5 coefficients of the 4-th order flux
            cf(0,i,im)=cosal(1)*eig11*skrt11*rkrt11+cosal(2)*eig21*skrt21*rkrt21 &
                      +eig11*cp011+eig21*cp021 
           cf(0,j,im)=cosar(1)*eig12*skrt12*rkrt12+cosar(2)*eig22*skrt22*rkrt22 &
                      +eig12*cp012+eig22*cp022
           cf(1,i,im)=3*(sinbl(1)*eig11*(ckrt11*rkrt11-skrt11*rksq11) &
                      +sinbl(2)*eig21*(ckrt21*rkrt21-skrt21*rksq21)) &
                      +eig11*cp111+eig21*cp121
           cf(1,j,im)=3*(sinbr(1)*eig12*(ckrt12*rkrt12-skrt12*rksq12) &
                      +sinbr(2)*eig22*(ckrt22*rkrt22-skrt22*rksq22)) &
                      +eig12*cp112+eig22*cp122
           cf(2,i,im)=5*(cosal(1)*eig11*rkrt11*(-3*ckrt11*rkrt11+3*skrt11*rksq11+skrt11) &
                      +cosal(2)*eig21*rkrt21*(-3*ckrt21*rkrt21+3*skrt21*rksq21+skrt21)) &
                      +eig11*cp211+eig21*cp221
           cf(2,j,im)=5*(cosar(1)*eig12*rkrt12*(-3*ckrt12*rkrt12+3*skrt12*rksq12+skrt12) &
                      +cosar(2)*eig22*rkrt22*(-3*ckrt22*rkrt22+3*skrt22*rksq22+skrt22)) &
                      +eig12*cp212+eig22*cp222
           cf(3,i,im)=7*(sinbl(1)*eig11*rksq11*(ckrt11*rkrt11*(15+ksq11)-15*skrt11*rksq11-6*skrt11) &
                      +sinbl(2)*eig21*rksq21*(ckrt21*rkrt21*(15+ksq21)-15*skrt21*rksq21-6*skrt21)) &
                      +eig11*cp311+eig21*cp321
           cf(3,j,im)=7*(sinbr(1)*eig12*rksq12*(ckrt12*rkrt12*(15+ksq12)-15*skrt12*rksq12-6*skrt12) &
                      +sinbr(2)*eig22*rksq22*(ckrt22*rkrt22*(15+ksq22)-15*skrt22*rksq22-6*skrt22)) &
                      +eig12*cp312+eig22*cp322
           cf(4,i,im)=9*(cosal(1)*eig11*rkrt11*(rksq11*(-5*ckrt11*rkrt11*(21+2*ksq11)+105*skrt11*rksq11 &
                      +45*skrt11)+skrt11)+cosal(2)*eig21*rkrt21*(rksq21*(-5*ckrt21*rkrt21*(21+2*ksq21) &
                      +105*skrt21*rksq21+45*skrt21)+skrt21))+eig11*cp411+eig21*cp421
           cf(4,j,im)=9*(cosar(1)*eig12*rkrt12*(rksq12*(-5*ckrt12*rkrt12*(21+2*ksq12)+105*skrt12*rksq12 &
                      +45*skrt12)+skrt12)+cosar(2)*eig22*rkrt22*(rksq22*(-5*ckrt22*rkrt22*(21+2*ksq22) &
                      +105*skrt22*rksq22+45*skrt22)+skrt22))+eig12*cp412+eig22*cp422
! Update the 4-th order source coefficients
            do l=0,4
              fsc(l,i)=fsc(l,i)+xsnf(i,im)*cf(l,i,im)
              fsc(l,j)=fsc(l,j)+xsnf(j,im)*cf(l,j,im)
            enddo
            enddo  ! for im

            do im=1,ng
              scs=0  
              do m=1,ng   
                do l=0,4
                  scs(l,i)=scs(l,i)+xssm(i,im)%from(m)*cf(l,i,m)      
                  scs(l,j)=scs(l,j)+xssm(j,im)%from(m)*cf(l,j,m)      
                enddo
              enddo  
              do l=0,4
                qs(l,i,im)=reigv*xchi(i,im)*fsc(l,i)+scs(l,i)          
                qs(l,j,im)=reigv*xchi(j,im)*fsc(l,j)+scs(l,j)             
              enddo
            enddo ! for im
! Current Convergence condition
          do im=1,ng
            j1_err(im)=abs(jnet(1,j,im)-j1_old(im))
            j2_err(im)=abs(jnet(2,j,im)-j2_old(im))
            j1_old(im)=jnet(1,j,im)
            j2_old(im)=jnet(2,j,im)
          enddo
          error1=maxval(j1_err)
          error2=maxval(j2_err)
          if (error1.lt.epsm6.and.error2.lt.epsm6) then
!           print*,"node=",i,j,"nodal=",iin
           exit
          endif 
! Tested which moment current is converged faster
! 2nd moment current converged fast in case of 1 box problem
! Random selection in case of multi-box problem
!            if (j11_err.lt.epsm5.and.j12_err.lt.epsm5) then
!              print*,iin,"0th moment"
!              exit
!            endif
!            if (j21_err.lt.epsm5.and.j22_err.lt.epsm5) then
!              print*,iin,"2nd moment"
!              exit
!            endif
          enddo ! for iin
!
          if (det.ne.0.and.j.eq.tnode) then
            call sanm1n(i,j)
          endif
        enddo ! for i
!      enddo ! for k

!#define cmfdp1
#ifdef cmfdp1
      do im=1,ng
        do i=1,tnode
          j=i+1          
          if (opt.eq.1) then
            avgflx(2,i,im)=-rxst(i,im)*rhz(i)*(0.4*(jnet(1,j,im)-jnet(1,i,im)) &
                          +0.6*(jnet(2,j,im)-jnet(2,i,im)))
            psphi(1,i,im)=avgflx(1,i,im)+2*avgflx(2,i,im)
          else
            exit
          endif
        enddo
!        avgflx(2,0,im)=avgflx(2,1,im)
!        avgflx(2,tnode+1,im)=avgflx(2,tnode,im)
      enddo
#endif

#ifdef m2
      do im=1,ng
        do i=1,tnode
          j=i+1
          jin(1,i,im)=0.25*sflux(1,i,im)+0.5*jnet(1,i,im)+0.3125*sflux(2,i,im)
          jin(2,i,im)=0.25*sflux(1,j,im)-0.5*jnet(1,j,im)+0.3125*sflux(2,j,im)
        enddo
      enddo
#endif

      end subroutine
!          jnet(1,i,im)=-beta11*((2+eig11)*krt11*(sinb(1,i,im)*ckrt11+cosa(1,i,im)*skrt11) &
!                       +(2+eig21)*krt21*(sinb(2,i,im)*ckrt21+cosa(2,i,im)*skrt21) &
!                       +eig11*(cp111+3*cp211+6*cp311+10*cp411)+eig21*(cp121+3*cp221+6*cp321+10*cp421) &
!                       +2*(cp111+cp121)+6*(cp211+cp221)+12*(cp311+cp321)+20*(cp411+cp421))           
!         jnet(2,i,im)=-beta21*(sinb(1,i,im)*ckrt11*krt11+sinb(2,i,im)*ckrt21*krt21 &
!                       +cosa(1,i,im)*skrt11*krt11+cosa(2,i,im)*skrt21*krt21+(cp111+cp121) &
!                       +3*(cp211+cp221)+6*(cp311+cp321)+10*(cp411+cp421))
!           low=1/((det1*det2)**2-(k2*k4-k1*k5)*(k11*k7-k10*k8) &
!             +det1*det2*(k10*k2-k11*k5-k1*k7+k4*k8))
!           up11=det2*(-k11*k3*k5+k11*k2*k6+k3*k4*k8-k1*k6*k8) &
!              +det1*det2*(-k12*k2+det2*k3+k1*k9)-(k2*k4-k1*k5)*(k12*k8-k11*k9)
!           up12=det2*(k10*k3*k5-k10*k2*k6-k3*k4*k7+k1*k6*k7) &
!              +(k2*k4-k1*k5)*(k12*k7-k10*k9)-det1*det2*(-k12*k5+det2*k6+k4*k9)
!           up21=-(k3*k5-k2*k6)*(k11*k7-k10*k8)+(det1**2)*det2*k9+det1*(-k12*k2*k7 &
!              +det2*k3*k7+k12*k5*k8-det2*k6*k8+k10*k2*k9-k11*k5*k9)
!           up22=-(det1**2)*det2*k12+(k3*k4-k1*k6)*(k11*k7-k10*k8) &
!              +det1*(-det2*k10*k3+det2*k11*k6+k1*k12*k7-k12*k4*k8-k1*k10*k9+k11*k4*k9)
